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VARIABLES 


Coordinates 


r 

plume,  radial  coordinate 

z 

plume  transverse  coordinate 

s 

primary  line-of-sight  coordinate 

0 

scattering  line-of-sight  coordinate 

0 

scattering  angle 

♦ 

azimuth  around  primary  line  of  sight 

Radial  Gas  Variables 

Pg(r) 

total  pressure 

cg(r) 

mole  fraction  concentration 

Vr> 

temperature 

k (r) 

band  model  absorption  parameter 

l/«(r) 

band  model  line  density  parameter 

Vr) 

pressure  broadened  (Lorentz)  line  halfwidth 

Vr) 

Doppler  line  halfwidth 

Bg(D 

Planck  function  at  Tg 

Radial  Particle  Variables 


°a(r) 

absorption  cross  section 

Os(r) 

scattering  cross  section 

dog(r,0)/dn 

differential  scattering  cross  section 

a(r) 

volume  absorption  cross  section 

B(r) 

volume  scattering  cross  section 

Y(r) 

volume  extinction  cross  section 
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p(r,e) 

scattering  phase  function 

Tp(r) 

temperature 

Bp(r) 

Planck  function  at  Tp 

Transverse  Emission/ Absorption/Scattering  Variables 


N(z) 

total  gaa-plus-particle  emission 

t(z) 

total  gas-plus-particle  transmission 

Vz) 

particle-only  emission 

Tp(z) 

particle-only  transmission 

f(0,z) 

laser  scattering  efficiency  function 

Calculation  Variables 

J(r) 

emission  source  function 

K(r) 

absorption  source  function 

G(z,r) 

emission  or  absorption  kernel  function 

T*(z,r) 

gas  transmittance  at  (s,r)  intersections 

y*(z,r) 

derivative  function  at  (s,r)  intersections 

Ta<z.r) 

transmittance  due  to  particle  absorption  at 

(s,r) 

intersections 

T*(z»r) 

transmittance  due  to  particle  scattering  at 

(s,r) 

intersections 

Qg(z.r) 

scattering  source  function  at  (s,r)  intersections 

J(r,0) 

laser  scattering  source  function 

G(z,r ,0) 

laser  scattering  kernel  function 

Miscellaneous 

X 

wavelength 

V 

wavenumber 

8 


plume  radius 

particle  radius 

particle  size  distribution 

complex  index  of  refraction 

real  part  of  index  of  refraction 

Imaginary  part  of  index  of  refraction 


1 .  INTRODUCTION 


ff 

EMABIC  is  a system  of  computer  codes  designed  for  infrared  radiation 
prediction  and  inversion  diagnostics  of  low-visibility  propellant  tactical 
rocket  motor  plumes.  The  prediction  capabilities  of  the  system  include  pro- 
files of  plume  infrared  emission  and  transmittance  transverse  to  the  plume 
axis  from  given  radial  profiles  of  the  plume  gas  and  particle  properties 
(e.g.,  concentration  and  temperature).  The  angular  distribution  of  laser 
radiation  scattered  by  the  particle  component  of  the  plume  is  also  pre- 
dicted. The  diagnostic  capability  is  the  retrieval  of  radial  gas  and  particle 
distributions  from  observed  transverse  and  laser  scattering  profiles. 

The  radiation  prediction  capabilities  of  the  system  derive  from  the 
previously  developed  prediction  code  EAPROF  (Refs.  1 and  2),  and  the  retrieval 
diagnostics  procedure  from  the  gas-only  plume  inversion  code  EMABIC  (Kef. 
3).  The  present  system  supersedes  both  of  these  codes.**  The  diagnostic 


EMABIC  is  an  acronym  for  Emission/Absorption  Inversion  J3odes. 

1.  S.  J.  Young,  Retrieval  of  Flow-Field  Gas  Temperature  and  Concentration 
of  Low-Visibility  Propellant  Rocket  Exhaust  Plumes,  AFRPL-TR-82-13 , 
U.  S.  Air  Force  Rocket  Propulsion  Laboratory,  Edwards  Air  Force  Base, 
California,  March  1982. 

2.  S.  J.  Young,  User's  Manual  for  the  Plume  Signature  Code  EaPROF,  AFRPL- 
TR-81-08,  U.  S.  Air  Force  Rocket  Propulsion  Laboratory,  Edwards  Air 
Force  Base,  California,  January  1981. 

3.  S.  J.  Young,  Inversion  of  Plume  Radiance  and  Absorptance  Data  for 
Temperature  and  Concentration,  AFRPL-TR-7M-6U,  6..  S.  Air  Force  Rocket 
Propulsion  Laboratory,  Edwards  Air  Force  Base,  California,  29  September 
1978. 

itft 

Except,  the  gas-only  diagnostics  code  with  error  propagation 

capability,  as  described  in  Ref.  4,  is  still  maintained. 

4.  S.  J.  Young,  Random  Error  Propagation  Analysis  in  the  Plume  Diagnostic 
Code  EMABIC,  AFRPL-TR-81-59 , U.  S.  Air  Force  Rocket  Propulsion 
Laboratory,  Edwards  Air  Force  Base,  California,  July  1981. 
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procedures  for  retrieving  particle  properties  and  an  improved  procedure  for 
gas  property  retrieval  with  account  of  particles  are  developed  in  Ref.  5. 
Reference  5 is  also  the  primary  companion  report  for  this  user’s  manual  and 
should  be  consulted  for  details  not  mentioned  here. 

The  primary  motive;'  on  for  EMABIC  was  the  development  of  tne  diagnos- 
tics tools  for  retrieving  the  radial  particle  and  gas  properties.  The  predic- 
tion capabilities  of  EMABIC  follow  as  procedures  required  in  the  iterative 
method  of  retrieval,  and  as  a means  of  generating  synthetic  data  with  which  to 
test  the  diagnostic  tools.  Application  of  the  diagnostic  scheme  requires 
traditional  transverse  infrared  emission/absorption  (E/A)  data  at  two  wave- 
lengths and  angular  scattering  data  at  one  of  the  wavelengths.  A^  is  a spec- 
tral position  at  which  both  the  gas  and  particle  species  of  interest  radiate 
and  absorb.  A^  is  a position  at  which  only  the  particles  radiate  and 
absorb.  The  transverse  E/A  data  at  these  wavelengths  are 


N(z)  i 

t(z)  ) 
Np(z)  | 

Tp(z)  j 


A = A 


1 


A * 


A 


2 


where  z(o  < z < k)  is  the  transverse  coordinate,  R is  the  plume  radius,  N is 
the  plume  radiance,  and  x is  the  transmittance  of  the  plume  to  external  radia- 
tion. Scattering  data  are  required  only  at  A^  and  consist  of  the  scattering 
intensity  efficiency  function 


b.  S.  J.  Young,  Retrieval  of  Flow-Field  Temperature  and  Concentration  ot 
Low-Visibility  Propellant  Rocket  Exhaust  Plumes,  AFKPL-TR-82-38,  U.  S. 
Air  Force  Rocket  Propulsion  Laboratory,  Edwards  Air  Force  Base, 
California,  February,  1983. 
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f(z,8) 


o < 9 < it 
o < z < R 


Thin  function  is  discussed  in  Kef.  5. 

From  these  data  profiles,  the  retrieval  codes  ofc  EMABIO  recover  the 
radial  plume  profiles  of 


a(r) 

particle  volume  absorption 

cross 

section 

8(r) 

particle  volume  scattering 

cross 

section 

Y(r) 

particle  volume  extinction 

cross 

section 

P(r,6) 

particle  scattering  phase  function 

T (r) 

P 

particle  temperature 

Vr> 

gas  temperature 

cg<r) 

gas  concentration 

The  prediction  codes  of  EMABIC  essentially  reverse  this  process  and  predict 
the  transverse  and  angular  scattering  data  from  the  radial  data. 

The  overall  structure  of  EMA&IC  is  described  in  Section  2.  Presently, 
the  system  consists  of  six  basic  driver  programs  and  45  subprogram  units.  The 
system  is  structured  for  easy  expansion  to  handle  envisioned  added  capability, 
some  of  which  are  mentioned  in  Section  2.  In  Section  3,  a detailed  descrip- 
tion or  data  preparation  for  the  six  basic  codes  Is  given.  An  example  appli- 
cation for  a minimum  smoke  propellant  plume  model  containing  1^0  and  AI2O3, 
respectively,  as  the  gas  and  particle  phases,  is  given  in  Section  4.  This 
example  is  taken  from  Ref.  5.  Most  of  the  computational  methods  employed  in 
the  system  have  been  discussed  in  the  references  to  its  predecessor  codes 
(Refs.  1,  2,  and  3).  One  significant  change  is  the  Abel  inversion  algo- 
rithm. The  revised  procedure  is  presented  in  the  Appendix. 
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2.  EMABIC  CODES  AND  STRUCTURE 


2.1  Overview 

EMABIC  was  conceived  as  a general  diagnostic  tool  for  the  retrieval  of 
gas  and  particle  properties  of  low-visibility  propellant,  tactical  rocket 
motor  plume0  from  experimental  E/A  and  scattering  measurements  made  on 
plumes.  The  program  name  DPREP1  has  been  reserved  for  a code  that  would 
process  such  data  (e.g.,  smoothing,  extrapolation,  scaling)  and  produce  the 
transverse  data  profiles  N(z),  t(z),  Np(z),  "^(z),  and  f(z,0)  required  by  the 
inversion  procedures.  In  lieu  of  experimental  data,  a set  of  transverse  data 
simulation  codes  has  been  envisioned  consisting  of  the  programs  DPR1.P2, 
EAPROF,  and  NOISE.  At  present,  only  the  programs  DPREP2  and  EAPROF  have  been 
completed.  DPREP2  is  a data  preparation  code  that  generates  a file  of  racial 
plume  properties  (e.g.,  temperature,  concentrations,  band  model  parameters). 
Its  primary  use  is  the  preparation  of  input  for  the  actual  transverse  profile 
generation  code  EAPROF.  NOISE  Is  a planned  program  that  will  add  random 
errors  to  the  profiles  generated  by  EAPROF  in  order  to  study  error  propagation 
In  the  retrieval  programs.  The  general  data  flow  for  the  preparation  of  input 
data  for  retrieval  is  shown  in  Fig.  1. 

The  inversion  programs  of  EMABIC  consist  of  three  programs  (PARIC1, 
PARIC2 , and  PARIC3)  for  the  retrieval  of  plume  particle  properties  anu  one 
code  (GASIC)  for  the  retrieval  of  plume  gas  properties.  The  general  flow  of 
data  for  retrieval  is  illustrated  in  Fig.  2.  Briefly,  PARIC1  retrieves  the 
volume  extinction  cross  section  y(r)  from  the  transverse  particle-only  trans- 
mittance profile  fp(z);  PARIC2  retrieves  the  volume  absorption  cross  section 
a(r),  the  volume  scattering  cross  section  3(r),  and  the  scattering  phase 
function  u(r,9)  from  the  laser  scattering  efficiency  function  f(z,0);  PAR1C3 
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retrieves  the  particle  temperature  profile  Tp(r)  from  the  particle-only  radi- 
ance profile  Np(z);  and  GASIC  retrieves  the  gas  temperature  Tg(r)  and  concen- 
tration c^(r)  from  the  coupled  gas-plus-particle  radiance  N(z)  and  transmit- 
tance t(z). 

A fourth  particle  retrieval  code  (PARIG4)  is  planned  that  would  allow 
the  retrieval  of  the  intrinsic  particle  scattering  cross  sections 
a (r),  a (r)  and  do  (r)/dfl  from  o(r),  0(r),  and  p(r,0),  respectively,  as 

9 8 S 

well  as  the  retrieval  of  the  absolute  particle  loading  np(r)  and  particle  size 
distribution  F(r,a)  (a  - particle  radius).  Retrieval  with  PAR1C1,  PAR1C2, 
PARIC3,  and  GASIC  requires  no  assumptions  on  the  index  of  refraction  or  shape 
of  the  plume  particles.  Implementation  of  PARIC4,  however,  would.  Its  devel- 
opment will  assume  that  the  particles  are  homogeneous  spheres  (so  that 
standard  Mie  scattering  theory  can  be  used)  and  that  the  particle  composition 
(and  thus  the  index  of  refraction  m ■ n -ix)  is  known. 

Data  flow  within  EMABIC  is  handled  with  three  random  access  files.  The 
file  generated  by  DPREP2  is  called  the  "a  priori  radial  data  file"  and  is  used 
primarily  as  input  for  the  transverse  profile  generation  program  EAPROF.  The 
output  of  EAPROF  is  written  to  a file  called  the  "transverse  data  file" 
(DPREP1  would  generate  this  same  file)  and  is  the  basic  input  for  the  re- 
trieval codes.  Finally,  the  retrieval  codes  write  their  results  to  a file 
called  the  "retrieved  radial  data  file.”  A summary  of  the  contents  of  these 
files,  and  the  codes  that  write  to  and  read  from  cheiu,  is  given  in  Table  1. 
Detailed  discussions  of  the  files  are  given  in  the  following  sections. 

In  general,  the  application  of  the  inversion  codes  to  actual  data  must 
be  made  in  the  order  PARIC1,  PARIC2,  PARIC3,  and  GASIC  because  the  retrieved 
result  of  each  code  is  required  by  the  subsequent  codes.  (PARIC4  could  follow 
Immediately  after  PARIC3;  the  results  of  GASIC  would  not  be  required.)  For 
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example,  Y(r),  which  is  obtained  with  PARIC1,  is  stored  on  the  "retrieved 
radial  data"  file  and  used  by  all  of  the  subsequent  codes.  Anticipating  that 
no  experiment  is  likely  to  be  performed  in  the  near  future  that  would  measure 
all  of  the  data  required  for  a full  implementation  of  the  EMABIC  retrieval 
capability,  provisions  have  been  made  for  operating  each  code  with  "assumed" 
results  from  a previous  code.  These  provisions  are  made  by  allowing  each  code 
to  obtain  its  required  radial  data  from  the  "a  priori  radial  data  file"  rather 
than  from  the  "retrieved  radial  data  file."  These  "assumed"  values  are  placed 
on  the  "a  priori  radial  data  file"  with  DPREP2. 

In  its  present  form,  EMABIC  consists  of  the  six  main  programs  DPREP2, 
EAPROF , PARICI,  PARIC2,  PARIC3,  and  GASIC,  45  subprogram  units,  and  three 
implied  data  files.  The  subprograms  form  a common  base  from  which  the  six 
main  programs  call  (either  directly  or  indirectly)  as  required.  Some  subpro- 
grams are  called  by  only  one  main  program,  but  most  are  called  by  more  than 
one.  Brief  descriptions  of  these  subprograms  are  given  in  the  following 
sections  as  they  arise  in  the  discussions  of  the  six  main  programs.  As  aids 
to  understanding  the  structure  of  EMABIC,  the  subprograms  are  grouped  accord- 
ing to  general  computational  function  in  Table  2,  and  the  subroutines  actually 
called  by  each  of  the  main  programs  are  indicated  in  the  routine-calling  flow 
diagrams  of  Figs.  3 through  8.  (These  should  not  be  interpreted  as  detailed 
calculation  flow  diagrams.) 

2.2  Radiation  Prediction  Codes 

The  new  codes  DPREP2  and  EAPROF  replace  the  older  code  EAPROF  (Refs.  I 
and  2)  and  add  the  capability  for  general  radial  variation  of  particle  scat- 
tering parameters  and  the  calculation  of  the  laser  scattering  efficiency 
function.  These  codes  compute  the  transverse  profiles  of  infrared  emission 
and  extinction  and  laser  scattering  for  an  axisymmetric,  axially  uniform, 
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Table  2.  Grouping  of  Subroutines  by  Computational  Function. 


Data  Preparation  Routines 

BPARAM  ZAFIT  NKPARAM* 

KDPARAM  ZONEFIT 

ANGLE IT 


Input/uutput  Routines  Convergence  Testing  Routines 


INPUT1 

OUTPUT 1 

C0NVRG1 

INPUT2 

OUTPUT2 

C0NVRG2 

INPUT3 

0UTPUT3 

C0NVRG3 

INPUT4 

OUTPUT4 

INPUTS 

OUTPUTS 

Radiation 

Routines 

SRCFUNC 

JSCAT 

PLANCK 

YCGL 

F 

KERNEL 

GSCAT 

YDRL 

YDRL 

G 

PROFILE 

FSCAT 

YLSL 

WMIX 

TRANS 

YCGD 

QSCAT 

YMLD 

RADNCE 

YLSD 

SLOS 

YMIX 

Inversion  Routines  Quadrature  Routines 


ABEL 

CTSOLVE 


SOLID 

QUAD 

COEFF 


Not  used  in  current  version  of  EMABIC. 

**These  ten  subroutines  comprise  the  subroutine  set  RADECK. 
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Fig.  3.  Subroutine-Calling  Flow  Diagram 
for  DPREP2.  Degree  of  Indenta- 
tion of  subroutine  name  indi- 
cates the  level  at  which  it  is 
called. 


22 


Fig.  4.  Subroutine-Calling  Flow  Dia- 
gram for  EAPROF.  Degree  of 
indentation  of  subroutine 
name  indicates  the  level  at 
which  it  is  called. 
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Fig.  5.  Subroutine-Calling  Flow  Diagram 
for  PARIC1.  Degree  of  indenta- 
| tion  of  subroutine  name  indi- 

cates the  level  at  which  it  is 
called. 
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cylindrical  plume  from  radial  profiles  of  gas  temperature,  pressure,  and 
concentration  and  particle  temperature  and  number  density.  The  radiation 
mc-iel  treats  gas  radiation  transfer  with  band  model  methods  and  particle 
radiation  transfer  with  the  single-scattering  approximation.  The  radiat  ton 
model  correctly  couples  the  gas  and  particle  components  into  a single 
emitting,  absorbing,  and  scattering  medium.  The  program  treats  just  one  gas 
species  and  one  particle  species  at  a time. 

The  gas  band  model  is  the  Malkmus  statistical  model  and  employs  either 
the  Curtis-Godson  (CG)  or  derivative  (DR)  approximations  to  handle  the  inhomo- 
geneity and  nonisothermality  of  the  plume.  Lorentz,  Doppler,  or  Voigt  gas 
spectral  line  profiles  may  be  used. 

The  single-scattering  geometry  used  for  particle  radiation  transport  is 
shown  in  Fig.  9.  The  s-axis  is  the  primary  line  of  sight  (LOS).  The  i,oS 
shown  in  Fig.  9 is  the  one  that  goes  through  the  full  plume  diameter.  Ac  ; tie 
LOS  is  scanned  out  across  the  lateral  extent  of  the  plume,  it  cuts  progres- 
sively shorter  chords  of  the  cylindrical  plume.  The  o-axis  is  the  scattering 
LOS.  It  is  described  oy  the  value  of  s where  it  branches  off  the  primary  LOS 
and  by  the  scattering  angle  9 and  the  azimuthal  angle  <f>.  The  single- 
scattering approximation  includes  radiation  emitted  along  the  primary  LOS  and 
radiation  that  has  been  scattered  once  from  the  scattering  LOS  into  the  pri- 
mary LOS.  If  the  scattering  LOS  terminates  on  the  nozzle  exit  plane,  motor 
radiation  scattered  into  the  primary  LOS  is  also  included.  The  exit  plane  is 
modeled  as  a solid  disc  with  uniform  temperature  and  emisslvity. 

Extinction  of  external  radiation  shone  through  the  plume  Is  assumed  to 
be  caused  by  gas  absorption,  particle  absorption,  and  particle  outscat- 
tering.  The  single-scattering  approximation  does  not  allow  inscattering  to 
contribute  to  extinction  calculations. 
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2.2.1  Program  DPKEF2 

The  purpose  of  DPKEP2  is  Co  prepare  a file  of  radial  plume  data. 
Various  data  defined  on  arbitrary  radial  and  angle  grids  are  read  in  by  DPREP2 
and  fitted  to  a desired  radial  grid,  scattering  angle  grid,  and  laser  scat- 
tering angle  grid.  The  radial  input  data  are  the  temperature  and  concen- 
tration profiles  for  the  gas  and  particle  species  of  interest,  the  band  model 

•r 

parameters  for  the  ga^r,  and  the  scattering  cross  sections  for  the  particles. 
The  input  grid  data  are  the  number  of  equai-size  radial  zones  to  which  the 
radial  data  are  fitted,  the  solid  angle  integration  scattering  angle  grid  0, 
and  the  laser  scattering  angle  grid  0 . Required  formats  for  these  data  are 
discussed  in  Section  3.2. 

The  first  step  of  data  preparation  is  the  fitting  of  the  input  tempera- 
ture and  concentration  data  for  both  the  gas  and  particles  to  the  equal-zone 
radial  grid.  This  is  performed  by  linear  interpolation  in  the  subroutine 
ZONEFIT.  With  the  gas  temperature,  pressure,  and  concentration  that  obtain  on 
the  grid,  quadratic  interpolation  is  made  on  the  input  band  model  parameters 
in  subroutine  KDPAKAM  in  order  to  obtain  the  parameters  on  the  radial  grid. 
Angle  independent  particle  scattering  parameters  are  fitted  to  the  radial  grid 
with  ZONEFIT  and  the  input  tabJe  of  scattering  parameters.  The  angle  depend- 
ent scattering  parameter  (the  differential  scattering  cross  section)  is  fitted 
to  the  radial  grid  and  the  integration  and  laser  scattering  angle  grids 

(0  and  0 )in  subroutine  ZAFIT.  This  subroutine  calls  ZONEFIT  for  the  re- 
L 

quired  radial  interpolation  and  ANGLFIT  tor  the  required  angle  interpola- 
tion. The  latter  interpolation  is  linear.  After  fitting  to  the  gr^ds 
0 and  9 , the  differential  scattering  cross  sections,  on  each  grid,  ire 

li 

renormalized  to  Air  by  use  of  subroutine  SOLID.  This  subroutine  calculates  the 
integral  over  4n  steradians  of  an  arbitrary  function  f(0)  defined  between 
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9*o  and  ir. 


After  further  straightforward  calculations,  the  a priori  radial  file  is 
written.  The  format  of  this  file  is  shown  in  Table  3. 

2.2.2  Program  EAPROF 

Program  EAPROF  calculates  transverse  E/A  and  laser  scattering  effici- 
ency functions.  Its  main  source  of  input  data  is  the  a priori  radial  data 
file.  These  data,  as  well  as  additional  data  in  card  form  (see  Section  3.3 
for  data  preparation)  are  read  by  1NPUT1.  The  basic  E/A  calculations  are 
performed  in  a loop  over  transverse  position  that  calls  the  three  subroutines 


SRCFUNC,  KERNEL,  and  PROFILE.  The  first  subroutine  computes  the  radial  source 

* 


function  J(r)  as  defined  in  Table  1 of  Ref.  5.  The  second  computes  the 
kernel  functions  G(z,r)  defined  in  the  same  reference.  The  third  subroutine 


computes  the  transverse  profiles  F(z)  defined  by  (see  Eq.  29,  Ref.  5) 


F(z)  -2  / J(r)  G(z,r)  rdr 


2 -2,1/2' 


, r - z 


The  numerical  quadrature  of  this  integral  is  performed  in  subroutine  ^UAD 
(with  quadrature  coefficients  computed  by  COEFF)  and  is  described  in  detail  in 
Ref.  3. 

Calculation  of  the  source  functions  J(r)  for  emission  requires  an 
evaluation  of  the  Planck  (blackbody)  radiation  function.  This  is  performed  in 
the  subprogram  PLANCK. 

Calculation  of  the  kernel  functions  G(z,r)  requires  the  prior  calcula- 
tion of  the  gas  transmittance  T*(z,r),  the  gas  derivative  functions  y^z.i), 


For  a full  appreciation  of  this  entire  section,  it  is  recommended  that 
Section  2 and  Table  1 of  Ref.  5 be  close  at  hand. 


Table  3.  A Priori  Radial  Data  File 


Index  Variable  Length 


Description 


1 

TITLE 

7 

DPREP2  job  run  title 

2 

GNAME 

1 

Gas  species  identification  name 

3 

PNAME 

1 

Particle  species  identification  name 

4 

AGRDID 

1 

Scattering  angle  grid  identification  name 

e 

LGRDID 

I 

Laser  scattering  grid  identification  name 

6 

GDTAID 

1 

Gas  data  identification  name 

7 

PDTAID 

1 

Particle  data  identification  name 

8 

GPRMID 

Gas  band  model  parameters  identification  name 

9 

PPRMID1 

1* 

Particle  index  of  refraction  identification  name 

!Q 

PPRMID2 

1 

Particle  scattering  parameters  identification  name 

ii 

V 

1 

1 

Wavenumber  of  gas  parameters  (cm  *) 

12 

V 

p 

1 

Wavenumber  of  particle  parameters  (cm  1 ) 

13 

Nr 

1 

Number  of  radial  grid  points 

14 

Ne 

1 

Number  of  scattering  angle  grid  points 

15 

V 

1 

Number  of  laser  scattering  angle  grid  points 

16 

L 

r 

51 

Radial  grid  points  (cm) 

17 

51 

Total  gas  pressure  (atm) 

18 

51 

Gas  temperature  (K) 

19 

TP 

51 

Particle  temperature  (K) 

20 

51 

Gas  concentration  (mole  fraction) 

21 

"p 

51 

Particle  concentration  (cm-J) 

22 

k 

51 

Absorption  band  model  parameter  (cm-1/atm) 

Line  density  band  model  parameter  (1/cm-1) 

23 

1/6 

51 

24 

yL 

51 

Pressure  broadened  line  width  (cm-1) 

25 

7d 

n 

51 

Doppler  broadened  line  width  (cm-1) 

26 

51 

Particle  index  of  refraction  (real  part) 

27 

K 

51* 

Particle  index  of  refraction  (imaginary  part) 
Particle  absorption  cross  section  (cnr) 

28 

a 

a 

51 

29 

0 

s 

51 

Particle  scattering  cross  section  (cm^) 

30 

a 

51 

Particle  volume  absorption  cross  section  (cm-1) 
Particle  volume  scattering  cross  section  (cm-1) 

31 

p 

51 

32 

9 

s 

37 

Scattering  angle  grid  (deg) 

33 

da  /dfi 
s 

1887 

Particle  differential  scattering  cross  section 
(cmVsr)  on  0_  grid 

Particle  scattering  phase  function  (sr-1)  on  0 

34 

Ps 

1887 

«r 

35 

< 

37 

Laser  scattering  angle  grid  (deg) 

36 

dOj  / di? 

1887 

Particle  differential  scattering  cross  section 
(cnr/sr)  on  P grid 

Particle  scattering  phase  function  (sr-1)  on  0^ 

37 

PL 

1887 

gr 

* 


Not  used  in  current  version 


of  EMABIC. 
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the  particle  transmittances  t *(z,r)  and  x^Cz.r),  and  the  scattering  source 

o p 

functions  Qg*(z,r).  All  except  Q*(z,r)  are  computed  in  subroutine  TRANS. 
This  subroutine  employs  the  band  model  procedures  described  in  Ref.  3 for 
computing  t (ztr)  and  y*(z,r).  In  addition  to  calls  to  QUAD  for  computing 
path  averaged  band  model  parameters,  any  of  the  secondary  radiation  routines 
of  the  set  KADECK  (i.e.,  YCGL,  YDRL,  YLSL,  YCGD,  YMLD,  YLSD,  YMIX,  F,  G,  or 
WMIXJ  may  be  called,  depending  on  the  gas  spectral  lineshape  and  path  nonuni- 
formity approximation  specified  as  part  of  the  EAPROF  input.  These  routines 
are  described  in  Ref.  3 and  works  referenced  therein.  The  particle  transmit- 
tance functions  T *(z,r)  and  T *(z,r)  are  calculated  directly  in  TRANS, 
a p 

The  scattering  source  functions  Qg*(z,r)  are  computed  in  subroutine 
QSGAT.  This  subroutine  calculates  the  integral  over  4x  steradians 
(o<9<n,  2ir)  of  radiation  impinging  on  a specified  point  on  an 

observation  line  of  sight.  The  integration  over  0 is  performed  in  subroutine 
SOLID.  The  actual  radiance  impinging  at  the  point  from  direction  (0,<J>)  is 
calculated  in  subroutine  RADNCE.  The  gas  and  particle  properties  along  the 
scattering  line  of  sight  required  for  the  radiance  calculation  are  set  up  in 
subroutine  SLOS.  The  gas  band  model  calculations  of  RADNCE  and  SLOS  parallel 
the  procedures  used  in  TRANS.  The  calculation  of  (^(z.r)  includes  the  con- 
tribution of  exit  plane  radiation. 

The  calculation  of  the  laser  scattering  efficiency  function  f(z,0) 
takes  place  in  a double  loop  over  laser  scattering  angles  0^  and  transverse 
position  z that  calls  the  subroutines  JSCAT,  GSCAT,  and  FSCAT.  The  first 
subroutine  computes  the  scattering  source  function  J(r,0)  (Eq.  43  of  Ret.  3), 
the  second  computes  the  kernal  function  G(z,r,0)  (Eq.  44  of  Ref.  5J,  and  the 
third  computes  the  scattering  efficiency  function  (Eq.  42  of  Ref.  5) 
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f(z,0)  - 2 / J(r,0)  G(z,r ,0)  — j~tr  . 

z (r  - z ) ' 2 

As  in  Che  subroutine  PROFILE,  this  integral  is  coapuCed  in  subroutine  QUAD. 

When  the  computation  of  the  E/A  and  scattering  functions  are  complete, 
the  results  are  listed  and  written  to  the  transverse  data  file  in  subroutine 
OUTPUT!.  The  format  of  the  output  file  is  given  in  Table  4. 

2. J Particle  Property  Retrieval  Codes 

2.3.1  Program  PARIC1 

Program  PARIC1  retrieves  the  radial  volume  extinction  profile  y(r)  from 
the  transverse  particle-only  transmittance  profile  t_(z).  The  transverse  data 
and  other  input  data  are  read  by  subroutine  INPUT2  (see  Section  3.4  for  data 
preparation).  Retrieval  is  accomplished  by  Abel  inversion  (Eq.  4d,  Kef.  5) 
of  fp(z)  fn  the  subroutine  ABEL  (see  Appendix).  The  result  for  y(r)  is 
written  to  the  retrieved  radial  data  file  by  0UTPUT2.  The  format  of  this  file 
is  given  in  Table  5. 

2.3.2  Program  PARIC2 

Program  PARIC2  retrieves  the  radial  volume  scattering  cross  section 
0(r),  volume  absorption  cross  section  a(r),  and  scattering  phase  function 
p(r,0)  from  the  transverse  laser  scattering  efficiency  function  f(z,0).  The 
transverse  data  and  other  input  data  are  read  in  by  subroutine  INPUT3  (see 
Section  3.5  for  data  preparation).  Retrieval  is  made  by  iterative  Abel  inver- 
sion for  Che  scattering  source  function  J(r,0)  (Eq.  49,  Ref.  5).  A first 

approximation,  J'(r,0),  is  made  for  J(r,6)  by  simple  Abel  inversion  of 

f(z,0).  This  first  approximation  is  then  used  to  compute  a new  function 
f'(z,0).  This  computation  is  made  exactly  as  is  done  in  EAPROF  (Section 

2.2.2)  except  chat  subroutine  JSCAT  is  not  called.  An  Abel  inversion  of  the 
difference  between  f and  f"  is  made  in  order  to  get  a correction  AJ  to 
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Table  4.  Transverse  Data  File 


Table  5.  Retrieved  Radial  Data  File 


- ' - f H 


Index 

Variable 

Length 

Description 

1 

TITLE 1 

7 

Transverse  data  file  TITLE 

2 

TITLE2 

7 

PARIC1  job  run  title 

3 

TITLE3 

7 

PARIC2  job  run  title 

4 

TITLE4 

7 

PARCI3  job  run  title 

5 

TITLE5 

7 

GASIC  job  run  title 

6 

TITLE6 

7* 

PARIC4  job  run  title 

7 

GNAME 

1 

Same  as  a priori  radial  data  file 

8 

PNAMK 

l 

Same  as  a priori  radial  data  file 

9 

Nr 

l 

Same  as  a priori  radial  data  file 

10 

Ke 

1 

Same  as  a priori  radial  data  file 

11 

NeL 

1 

Same  as  a priori  radial  data  file 

12 

No 

1* 

Number  of  particle  sizes 

13 

r 

51 

Same  as  a priori  radial  data  file 

14 

Tg 

51 

Same  as  a priori  radial  data  file 

15 

rl 

51 

Same  as  a priori  radial  data  file 

16 

51. 

Same  as  a priori  radial  data  file 

17 

np 

51 

Same  as  a priori  radial  data  file 

18 

a 

51 

Same  as  a priori  radial  data  file 

19 

p 

51 

Same  as  a priori  radial  data  file 

20 

Y 

51. 

Volume  extinction  cross  section  (cm"*) 

21 

a 

a 

51 

Same  as  a priori  radial  data  file 

22 

0 

s 

51* 

Same  as  a priori  radial  data  file 

23 

e 

s 

37 

Same  as  a priori  radial  data  file 

24 

da  /(W 
s 

1887* 

Same  as  a priori  radial  data  file 

25 

p8 

1887 

Same  as  a priori  radial  data  file 

26 

6L 

37 

Same  as  a priori  radial  data  file 

27 

da^/ d0 

1887* 

Same  as  a priori  radial  data  file 

?n 

1*1. 

1887 

Same  as  a priori  radial  data  file 

29 

a 

30* 

Particle  size  grid  (pm) 

30 

F 

1530* 

Particle  distribution  (pm  ) 

J'(r,e).  If  AJ  is  small  enough  (as  determined  by  an  input  convergence  cri- 
terion and  subroutine  C0NVRG1),  iteration  ceases.  Otherwise,  it  is  continued 
with  J'  replaced  by  J"  + AJ.  When  J(r,0)  is  determined  to  within  the  conver- 
gence criterion,  it  is  Integrated  over  4ir  steredlans  in  subroutine  SOLID  to 
obtain  0(r)  (Eq.  50,  Ref.  5).  Then,  a(r)  » y(r)  0(r)  and  p(r,0)  * 

4x  J(r,0)/g(r).  The  results  are  written  to  the  retrieved  radial  data  file  in 
OUTPUT3 . 

2.3.3  Program  PARIC3 

Program  PARIC3  retrieves  the  radial  particle  temperature  profile 
Tp(r)  from  the  transverse  particle-only  radiance  Np(z).  The  transverse  and 
other  input  data  are  read  in  by  subroutine  INPUT4  (see  Section  3.6  for  data 
preparation).  In  addition  to  reading  the  data,  INPUT4  also  fits  and  renormal- 
izes the  scattering  phase  function  obtained  in  PARIC3  (which  is  defined  on  the 
laser  scattering  grid  0L)  to  an  integration  scattering  angle  grid  0.  The 
fitting  and  renormalization  are  performed  with  the  subroutines  ANGLFIT  and 
SOLID.  Tha  actual  retrieval  of  Tp(r)  occurs  by  iterative  Abel  inversion  of 

Np(z)  to  obtain  the  source  function  J(r)  ■ a(r)  Bp(r)  (Eq.  55,  Ref.  5).  The 
iterative  computation  of  Np(z)  is  performed  as  in  EAPROF  except  that  the 
source  function  subroutine  SRCFUNC  is  not  used.  Convergence  of  the  iteration 

Is  tested  Jn  CONVRG2.'  When  J(r)  (and  thus  B_(r)J  is  determined  to  within  the 

F 

convergence  criterion,  a straightforward  inversion  of  the  Planck  function  is 
made  for  TI}(r).  The  result  is  written  to  tnt  ’•etrieved  radial  data  file  in 
WPUT4 . 

2.4  Gas  Property  Retrieval  Code  (Program  GASIC) 

Program  GASIC  retrieves  the  radial  gas  temperature  T (r)  and  concentra- 

o 

ticm  v-  g(r)  from  the  transverse  total  radiance  N(z)  and  transmittance  t(z) 
profiles.  The  transverse  and  other  data  are  read  in  by  subroutine  INPUT5  («fe 
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Section  3.7  for  data  preparation).  Retrieval  is  made  by  simultaneous  iter- 
ative Abel  inversion  of  the  total  transverse  radiance  N(z)  and  transmit- 


tance t(z)  profiles  to  obtain  the  source  functions  ic(r)  B (r)  and  <(r)  (E<j. 

8 

5b,  Ref.  5),  At  each  iteration,  the  subroutine  CTSOLVE  is  used  to  retrieve 
T (r)  and  c2(r)  from  the  current  source  functions,  and  convergence  is  tested 

o o 

in  subroutine  C0NVRG3.  The  iterative  computation  of  N(z)  and  t(z)  is  per- 
formed as  in  EAPROF  except  that  the  source  function  subroutine  SRCFUNC  is  not 
called.  With  each  iteration  result  for  Tff  and  c„,  new  values  for  band  model 

o o 

profiles  on  the  radial  grid  are  determined  in  subroutine  BPARAM  (which  calls 

KDPARAM).  When  Tg  and  Cg  are  determined  to  within  the  convergence  criteria, 

the  results  are  written  to  the  retrieved  radial  data  file  in  0UTPUT5. 

The  inversion  procedure  described  above  is  for  the  full  retrieval 

scheme  developed  in  Ref.  5.  Slight  modifications  are  made  for  the  gas-only 

and  first-order,  off-band  retrieval  approximations  options.  In  the  first 

approximation,  a gas-only  inversion  is  made  on  N(z)  and  t(z)  as  though  they 

were  for  a purely  gaseous  plume.  Particle  effects  are  totally  neglected,  in 

the  second,  a gas-only  inversion  is  made  on  N(z)  - N (z)  and  t(z)/t  (z). 

P P 


I: 
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3.  PREPARATION  OF  INPUT  DATA 


3.1  General  Instructions 

A computational  run  of  each  of  the  programs  DPREP2,  EAPROF,  PARIC1 , 
PARIC2,  PARIC3 , OR  GASIC  requires  a set  of  program  control  cards  to  specify 
the  mode  of  computation  and  to  supply  input  data.  Some  program  control  caras 
simply  specify  a computation  mode,  some  specify  a computation  mode  and  supply 
data,  while  others  signal  the  codes  that  blocVs  of  auxiliary  data  are  now  to 
be  read.  Each  type  of  control  card  contains  an  alphanuoeric  name  in  the  first 
ten  card  columns.  These  names  must  be  spelled  correctly  and  must  be  lei't- 
justiflei.  if  data  are  specified  on  a program  control  card,  they  must  be 
entered  in  accordance  with  the  format  specification  indicated  in  the  detailed 
description  of  each  card  given  below.  All  fields  of  the  program  control  cards 
are  10  columns  wide.  In  general,  integer  and  alphanumeric  data  must  be  right- 
ju8tifiad  in  their  fields.  Noninteger  numerical  data  may  be  entered  in  either 
F or  E formats  (with  decimal  point  and,  for  the  latter,  the  exponential  symbol 
E).  E-formatted  data  must  be  right-justified  in  their  fields.  These  same 
rules  aaply  to  data  entered  on  auA*  liary  card  decks. 

lach  of  the  programs  has  a multiple  run  feature.  After  all  of  the  data 
required  for  a run  have  been  entered  by  way  of  control  earns,  ayvniary  card 
decks,  cr  attached  data  files,  execution  is  initiated  with  the  RUN  control 
card.  When  the  computation  and  results  listing  are  completed,  the  program 
continues  to  read  program  control  cards  until  a new  RUN  card  is  encountered. 
A new  computation  is  then  begun  for  all  of  the  conditions  and  data  of  the 
first  run  except  those  that  have  been  changed  by  the  intervening  program 
control  cards  and  auxiliary  data  decks.  This  process  is  repeated  until  an 
end-of-file  card  is  encountered.  With  this  feature,  a large  number  of  related 
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runs  can  be  made  with  one  job  submission.  All  of  the  codes  write  their  re- 
sults to  a data  file  (TAPE2,  TAPE3,  OR  TAPE4).  If  multiple  runs  are  made  for 
any  of  the  codes  with  one  job  submission,  only  the  results  of  the  last  run  are 
saved  on  the  appropriate  file. 

In  general,  all  program  control  cards  that  call  for  the  read  in  of  data 
from  card  decks  or  attached  files  have  a print  option  for  listing  the  data. 
If  the  variable  PRINT  (format  A10)  on  these  cards  has  the  alphanumeric  value 
PRINT,  the  data  will  be  listed.  The  description  of  this  option  is  not  re- 
peated in  the  following  sections. 

A control  card  common  to  all  of  the  program  is  the  title  card.  Its 
name  is  TITLE,  and  columns  II  through  80  of  this  csrd  may  be  used  for  any  run 
identification  title  desired.  Use  of  the  card  is  optional.  Again,  the 
description  of  this  card  Is  not  repeated  in  the  following  sections. 

Other  than  the  requirement  that  all  required  data  be  specified  before  a 
RUN  card  Is  encountered  and  that  auxiliary  data  on  card  decks  immediately 
follow  the  control  card  that  calls  for  them,  the  program  control  cards  may  for 
the  most  part,  be  arranged  in  any  order.  Exceptions  are  specifically  men- 
tioned in  the  detailed  descriptions  that  follow. 

Great  care  should  be  taken  in  the  preparation  of  input  data  since  very 
few  checks  of  data  consistency,  setting  of  default  values,  or  issuance  of 
error  messages  are  provided.  A general  feature  of  data  preparation  is  that, 
if  particular  data  on  a card  are  not  required,  they  need  not  be  specified.  If 
none  cf  the  data  on  a control  card  is  needed,  that  card  need  not  be  included. 
3.2  DPREP2  Data  Preparation 

Program  DPREP2  prepares  radial  and  angular  data  for  use  by  the  E/n 
prediction  code  EAPROF.  The  types  of  control  cards  required  for  a run  oi 
program  DPREP2  are  illustrated  in  Fig.  10.  A description  of  each  type 


follows. 
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Title  Card 


f • 


t 


1. 

2.  Miscellaneous  Data  Card.  The  card  name  is  DATA.  NZONES  (format 
110)  is  the  number  of  equal-size  radial  zones  into  which  the 
plume  is  divided.  The  maximum  allowed  value  is  50.  (NZONES 
also  defines  the  number  of  transverse  zones  used  in  EAPROF.) 
GCOL  and  POOL  (each  format  110)  identify,  respectively,  which 
gas  and  particle  species  are  to  be  selected  from  the  auxiliary 
input  decks  of  radial  gas  and  particle  data  (Figs.  12  and  13, 
respectively).  GCOL  * 1,  2,  3,  or  4;  PCOL  * 1,  2,  or  3. 

3.  Scattering  Angle  Integration  Grid  Card.  The  card  name  is 
AGR1D.  This  card  calls  for  the  read  in  of  the  scattering  angle 
grid  that  EAPROF  will  use  for  Integration  over  solid  angles. 
The  required  deck  structure  for  these  data  is  shown  in  Fig.  11. 

4.  Laser  Scattering  Angle  Grid  Card.  The  card  name  is  LGRXD.  This 
card  calls  for  the  read  in  of  the  angle  grid  for  which  EAPROF 
will  compute  the  laser  scattering  efficiency  function  f(z,0). 
The  required  deck  structure  for  the  data  is  shown  in  Fig.  11. 

5.  Radial  Gas  Data  Card.  The  card  name  is  GDATA.  This  card  calls 
for  the  read  in  of  radial  gas  data.  The  required  deck  structure 
for  the  data  is  shown  in  Fig.  12. 

6.  Radial  Particle  Data  Card.  The  card  name  is  PDATA.  This  card 
calls  for  the  read  in  of  radial  particle  data.  The  required 
deck  structure  for  the  data  is  shown  in  Fig.  13. 

7.  Gas  Band  Model  Parameters  Card.  The  card  name  is  GPARAM.  This 
card  cal  Is  for  the  read  in  of  band  model  parameters  for  the  gas 
species  of  interest.  It  is  the  user's  responsibility  to  ensure 
that  the  parameters  read  in  are  consistent  with  the  gas  species 
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asked  for  by  GCOL  on  card  type  2.  The  required  deck  structure 
for  the  data  is  shown  in  Fig.  14. 

8,  Particle  Parameters  Card.  The  card  name  is  PPARAM2 . * This  card 
calls  for  the  read  in  of  particle  scattering  parameters  for  the 
particle  species  of  interest.  It  is  the  user’s  responsibility 
to  ensure  that  the  parameters  read  in  are  consistent  with  the 
particle  species  asked  for  by  PCOL  on  card  type  2.  The  required 
deck  structure  for  the  data  is  defined  by  the  FORTRAN  read 
routine  shown  in  Fig.  15. 

9.  Data  Listing  Card.  The  card  name  is  LIST.  The  default  output 
of  DPREP2  is  the  a priori  radial  data  file  TAPE2.  A listing  of 
the  prepared  data  is  made  if  the  control  card  LIST  is  included 
in  the  control  card  deck. 

10.  Listing  Suppression  Card.  The  card  name  is  NOLIST.  Listings 
for  subsequent  runs  can  be  suppressed  with  the  NOLIST  card  if 
listing  was  enabled  in  previous  runs  with  a LIST  card. 

11.  Execution  Card.  The  card  name  is  RUN.  When  this  card  is  en- 
countered, computations  are  begun  using  the  data  entered  up  to 
that  point,  an  output  listing  of  the  results  is  made,  and  the  a_ 
priori  radial  data  file  TAPE2  is  generated. 

3.3  EAPROF  Data  Preparation 

Program  EAPROF  calculates  transverse  E/A  and  laser  scattering  effici- 
ency functions.  Its  main  source  of  input  data  is  the  a priori  radial  data 
file  TAPE2  previously  generated  by  DPREP2.  Additional  data  required  to  run 


Card  name  PPAKAHl  is  reserved  for  future  use  to  read  in  particle  index 
of  refraction  data. 
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READ  (5,  100)  PPRMID,  NR,  NA,  WN 
READ  (5,  101)  (R(I),  I « 1,  NR) 
READ  (5,  101)  (SA(I),  I « 1,  NR) 
READ  (5,  101)  (SS(I),  1-1,  NR) 
READ  (5,  101)  (ANG(J),  J =*  1,  NA) 
DO  1 I = 1,  NR 

READ  (5,  101)  (DS(I,J) , J - 1,  NA) 
1 CONTINUE 

100  FORMAT  (A10,  2110,  E10.0) 

101  FORMAT  (8E10.0) 


PPRMID 

NR 

NA 

WN 

R(I) 

SAU) 

SS(I) 

ANG(J) 


Particle  parameters  identification  name. 

NumDer  of  radial  points  (NR  < 201). 

Number  of  angles  (NA  < 181). 

Wavenumber  ( cm~^ ) . 

Radial  positions  r(cm).  0 » R(l)  < ...  < r(nr)  = r. 
Absorption  cross  section  a (cm2)  at  R(I). 

Total  scattering  cross  section  a (cm2)  at  R(I). 

Scattering  angles  6 (deg).  0 - ANG(1  < ...  < ANG(NA)  » 18U. 


DS(1,J)  ^g[jJential  scatterin8  cross  section  do^dfi  (cm2/sr)  at  R(I)  and 


Fig,  15.  Input  Routine  for  Particle  Scattering  Parameters 
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EAPROF  are  supplied  by  the  control  cards  shown  in  Fig.  16  and  described  below. 

1.  Title  Card. 

2.  Calculation  Data  Card.  The  card  name  is  CALCDATA.  The  variable 
MODE  (format  A10)  determines  the  gas  and  particle  case  for  which 
calculations  will  be  made.  If  MODE  has  the  value  GAS,  calcula- 
tions will  be  made  as  though  only  the  gas  species  were  optically 
active.  The  value  PARTICLE  implies  that  only  the  particle 
species  is  active.  The  value  BOTH  implies  a coupled  gas-plus- 
particle  calculation.  If  MODE  has  the  value  ALL,  calculations 
are  performed  for  all  three  cases.  The  variables  SHAPE  and 
INHOM  (each  format  A10)  specify,  respectively,  the  gas  absorp- 
tion lineshape  and  optical  path  nonuniformity  approximation  to 
be  employed  in  the  gas  band  model.  SHAPE  must  be  one  of  the 
values  LORENTZ,  DOPPLER,  or  VOIGT.  INHOM  must  have  either  the 
value  CG  (for  the  Curtis-Godson  approximation)  or  DR  (for  the 
derivation  approximation).  NSCAT  (format  110)  is  the  number  of 
equal-length  segments  that  a scattering  LOS  is  divided  into  for 
numerical  integration.  Its  maximum  allowed  value  is  100.  NAZ 
(format  110)  is  the  number  of  intervals  over  which  the  360° 
azimuthal  angle  integration  is  performed  in  integrations  over 
solid  angle.  There  is  no  dimensional  limit  to  the  value  of  NAZ. 

3.  Plume  Data  Card.  The  card  name  is  PLMDATA,  LI  (format  E10)  is 
the  distance  (cm)  from  che  nozzle  exit  plane  to  the  observation 
scanning  plane.  L2  (format  E10)  Is  the  distance  (cm)  from  the 
scanning  plane  to  the  end  of  the  plume.  TN  (format  E10)  is  the 
temperature  (K)  of  the  nozzle  exit  plane  disc,  and  EN  (format 
E10)  is  its  emissivity.  These  data  (and  thus  the  card)  are  not 
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required  if  MODE  = GAS  on  the  calculation  data  card  (card  type 

2). 

4.  Radial  Data  Card.  The  card  name  is  RDATA.  This  card  calls  for 
the  read  in  of  the  a priori  radial  data  file.  This  file  must  be 
attached  to  the  EAPROF  job  run  as  TAPE2. 

5.  Execution  Card.  The  card  name  is  RUN.  When  this  card  is  en- 
countered, computations  are  begun  using  the  data  entered  up  to 
that  point,  an  output  listing  of  the  results  is  made,  and  the 
transverse  data  file  TAPE3  is  generated. 

3.4  PAKIC1  DATA  PREPARATION 

Program  PARIC1  retrieves  the  radial  volume  extinction  profile  y(r)  by 
Abel  inversion  from  tne  transverse  particle-only  transmittance  t^(z).  This 
transmittance  is  read  in  from  the  transverse  data  file  TAPE3  generated  by 
EAPROF  or  DPREP1.  The  control  cards  for  PARIC1  are  shown  in  Fig.  17  and 
described  below. 

1 . Title  Card. 

2.  Transverse  Data  Card.  The  card  name  is  ZDATA.  This  card  calls 

for  the  read  in  of  t (z)  from  the  transverse  data  file.  This 

P 

file  must  be  attached  to  the  PARIC1  job  run  as  TAPE3. 

3.  Execution  Card.  The  card  name  is  RUN.  When  this  card  is  en- 
countered, computations  are  begun  using  the  data  entered  up  to 
that  point,  an  output  listing  of  the  results  is  made,  and  the 
retrieved  function  y(r)  is  written  to  the  retrieved  radial  data 


file  TAPE4 


TITLE 

ZDATA 

RUN 
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PAK1C2  Data  Preparation 

Program  PARIC2  retrieves  the  radial  volume  scattering  cross  section 
$(r),  volume  absorption  cross  section  a(r),  and  scattering  phase  function 
p(r,0)  by  iterative  Abel  inversion  from  the  transverse  laser  scattering 
efficiency  function  f(z,0).  This  function  is  read  in  from  the  transverse  data 
file  TAPE3  generated  by  EAPROF  or  DPREPl.  The  control  cards  for  PAKIC2  are 
shown  in  Fig.  18  and  described  below. 

1 . Title  Card. 

2.  Convergence  Data  Card.  The  card  name  is  CONVERGE.  This  card 
enters  the  iterative  Abel  inversion  convergence  criteria.  IMaX 
(format  110)  is  the  maximum  number  of  iterations  allowed  (un- 
limited), and  ERROR  (format  E10)  is  the  maximum  rms  difference 
(percent)  allowed  between  successive  iterations  of  the  retrieved 
scattering  source  function  for  convergence  to  be  deemed 
complete. 

3.  Miscellaneous  Data  Card.  The  card  name  Is  NDATA.  NS CAT  (format 
110)  is  the  number  of  equal-length  segments  that  a scattering 
LOS  is  divided  into  for  numerical  integration.  Its  maximum 
allowed  value  is  100.  If  NL  (format  110)  has  the  value  1, 
calculations  will  be  performed  as  though  the  laser  scattering 
efficiency  function  were  defined  only  for  a single  line  of  sight 
through  the  plume  at  z = o.  The  retrieved  phase  function  and 
volume  scattering  cross  section  will  be  constant  in  r.  If 
NL  * 1 , a full  inversion  for  the  radial  variation  of  these  func- 
tions is  made. 

4.  Iteration  Listing  Card.  The  card  name  is  LIST.  Normally,  the 
intermediate  results  for  current  phase  function  and  transverse 


profiles  generated  at  each  iteration  are  not  listed.  The  LIST 
control  card  can  be  used  to  obtain  them. 

5.  Listing  Suppression  Card.  The  card  name  is  NOLIST.  Listing  ot 
intermediate  iteration  results  can  be  suppressed  with  the  NOLIST 
card  if  listing  was  enabled  in  previous  runs  with  a LIST  card 
(card  type  4). 

6.  Transverse  Data  Card.  The  card  name  is  ZDATA.  This  card  calls 
for  the  read  in  of  f(z,9)  from  the  transverse  data  file.  This 
file  must  be  attached  to  the  PARIC2  job  run  as  TAPE3. 

7.  Radial  Data  Card.  The  card  name  is  KDATA.  This  card  calls  for 
the  read  in  of  radial  data.  The  only  radial  data  required  by 
this  code  is  the  volume  extinction  cross  section  y(r).  If  CFILE 
(format  110)  has  the  value  2,  y(r)  is  obtained  from  the  a priori 
radial  data  file  TAPE2.  If  GFILE  has  the  value  4,  y(r)  is 
obtained  from  the  retrieved  radial  data  file  TAPE4.  The  appro- 
priate file  must  be  attached  to  the  PARIC2  job  run. 

8.  Execution  Card.  The  card  name  is  RUN.  When  this  card  is  en- 
countered, computations  are  begun  using  the  data  entered  up  to 
that  point,  an  output  listing  of  the  results  is  made,  and  the 
retrieved  functions  o(r),  0(r)  and  p(r,0)  are  written  to  the 
retrieved  radial  data  file  TAPE4,  p(r,9)  is  defined  on  the 
laser  scattering  angle  grid  9 . 

Lt 

3.6  PARIC3  Data  Preparation 

Program  PARIC3  retrieves  the  radial  particle  temperature  profile  T^lr) 
by  iterative  Abel  inversion  from  the  transverse  particle-only  radiance  profile 
Np(z).  This  radiance  is  read  in  from  the  transverse  data  file  TAPE3  generated 
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by  EAPROF  or  DPRKP1.  The  control  card  formats  for  PARIC3  are  shown  In  Fig.  19 
and  described  below* 

1.  Title  Card. 

2.  Convergence  Data  Card.  The  card  name  is  CONVERGE.  This  card 
enters  the  iterative  Abel  inversion  convergence  criteria.  IMAX 
(format  110)  is  the  maximum  number  of  iterations  allowed  (un- 
limited), and  ERROR  (format  E10)  is  the  maximum  rms  difference 
(percent)  allowed  between  successive  iterations  of  the  retrieved 
thermal  source  function  for  convergence  to  be  deemed  complete. 

3.  Plume  Data  Card.  The  card  name  is  PLMDATA.  LI  (format  E10)  is 
the  distance  (cm)  from  the.  nozzle  exit  plane  to  the  observation 
scanning  plane.  L2  (format  E10)  is  the  distance  (cm)  from  the 
scanning  plane  to  the  end  of  the  plume.  TN  (format  E10)  is  the 
temperature  (K)  of  the  nozzle  exit  plane  disc,  and  EN  (format 
E10)  is  its  emissivlty. 

A.  Miscellaneous  Data  Card.  The  card  name  is  NDATA.  NSCAT  (format 
110)  is  the  number  of  equal-length  segments  that  a scattering 
line  of  sight  is  divided  into  for  numerical  integration.  Its 
maximum  allowed  value  is  100.  NAZ  (format  110)  is  the  number  ol 
intervals  over  which  the  360°  azimuthal  angle  integration  is 
performed  in  integrations  over  solid  angle.  There  is  no  dimen- 
sional limit  to  NAZ. 

5.  Iteration  Listing  Card.  The  card  names  is  LIST.  Normally,  the 
Intermediate  results  for  current  thermal  source  function  and 
transverse  profiles  generated  at  each  iteration  are  not 


listed.  The  LIST  control  card  can  be  used  to  obtain  them 


6.  Listing  Suppression  Card.  The  card  name  is  NOLIST.  Listing  of 
intermediate  iteration  results  can  be  suppressed  with  the  NOLIST 
card  if  listing  was  enabled  in  previous  runs  with  a LIST  card 
(card  type  5). 

7.  Scattering  Angle  Integration  Grid  Card.  (Type  1)  The  card  name 
is  AGRID.  If  this  card  is  used,  it  must  come  before  the  Ri)ATA 
control  card  (card  type  10)  in  the  control  card  deck.  The  AGRID 
card  calls  for  the  read  in  of  the  scattering  angle  grid  0 that 
PARIC3  will  use  for  integration  over  solid  angle.  The  required 
'.deck  structure  for  these  data  is  shown  in  Fig.  11.  If  no  AGK1L) 
card  is  used,  the  scattering  grid  0 for  solid  angle  integration 

V 

and  the  values  of  the  scattering  phase  function  p on  that  grid 
will  be  read  (when  the  RDATA  card  is  encountered)  from  the  ji 
priori  radial  data  file  TAPE2.  If  an  AGRID  card  is  used,  the 
laser  scattering  grid  0 and  the  values  of  the  scattering  phase 
function  p^  on  that  grid  will  be  read  (when  the  RDATA  card  is 
encountered)  from  the  a^  priori  radial  data  file  TAPE2  or  the 
retrieved  radial  data  file  TAPE4,  and  an  interpolation  made  to 
get  the  scattering  phase  function  p on  the  grid  0 read  in  by  the 
AGRID  card. 

8.  Scattering  Angle  Integration  Grid  Card  (Type  2).  The  card  name 
is  NOAGRID.  This  control  card  restores  the  scattering  grid 
condition  of  subsequent  runs  to  the  equivalent  state  of  never 
having  encountered  an  AGRID  card. 

9.  Transverse  Data  Card.  The  card  name  is  ZDATA.  This  card  calls 
for  the  read  in  of  Np(z)  from  the  transverse  data  file.  This 
file  must  be  attached  to  the  PAKIC3  job  run  as  TAPE3. 
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10.  Radial  Data  Card.  The  card  name  is  RDATA.  This  card  calls  for 


the  read  In  of  radial  data.  The  radial  data  required  by  PARIC3 
are  the  volume  absorption  and  scattering  cross  sections, 
a(r)  and  $(r),  respectively,  and  the  scattering  phase  function 
p(r,0).  If  AFILE  (format  110)  has  the  value  2,  the  a(r)  data 
is  obtained  from  the  a priori  radial  data  file  TAPE2.  If  AFILE 
has  the  value  4,  a(r)  is  obtained  from  the  retrieved  radial  data 
file  TAPE4.  BFILE  and  PFILE  control,  respectively,  the  input 
of  B(r)  and  p(r,0)  in  the  same  way.  The  appropriate  files  must 
be  attached  to  the  job. 

11.  Execution  Card.  The  card  name  is  RUN.  When  this  card  is  en- 
countered, computations  are  begun  using  the  data  entered  up  to 
that  point,  an  output  listing  of  results  is  made,  and  the  re- 
trieved function  Tp(r)  is  written  to  the  retrieved  radial  data 
file  TAPE4. 

3.7  GASIC  Data  Preparation 

Program  GASIC  retrieves  the  radial  gas  temperature  T (r)  and  concentra- 

O 

tion  Cg(r)  by  coupled  iterative  Abel  inversion  from  the  transverse  total 
radiance  N(z)  and  transmittance  x(z)  profiles.  These  transverse  functions  are 
read  in  from  the  transverse  data  file  TAPE3  generated  by  EAPROF  or  DPREP1. 
The  control  card  formats  for  GASIC  are  shown  in  Fig.  20  and  described  below. 

1.  Title  Card. 

2.  Calculation  Data  Card.  The  card  name  is  CALCDATA.  The  variable 
MODE  (format  A10)  determines  the  inversion  mode.  If  MODE  has 
the  value  GAS,  an  inversion  is  made  as  though  the  input  trans- 
verse profiles  N(z)  and  t(z)  were  for  a purely  gaseous  plume. 
If  MODE  has  the  value  F00B,  the  first-order,  off-band  correction 
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procedure  is  used  for  inversion , and  if  MODE  has  the  value  BOTH, 
the  fully-coupled,  inversion  procedure  is  used*  The  variables 
SHAPE  and  1NH0M  (each  format  A10)  specify,  respectively,  the  gas 
absorption  lineshape  and  optical  path  nonuniformity  approxima- 
tion to  be  employed  in  the  gas  band  model.  SHAPE  must  be  one  of 
the  values  LORENTZ,  DOPPLER,  or  VOIGT.  INHOM  must  have  either 
the  value  CG  (for  the  Curtis-Godson  approximation)  or  DR  (for 
the  derivative  approximation). 

3.  Convergence  Data  Card.  The  card  name  is  CONVERGE.  This  card 
enters  the  iterative  Abel  inversion  convergence  criteria.  IMAX 
(format  110)  is  the  maximum  number  of  iterations  allowed  (un- 
limited). TERROR  and  CERROR  (each  format  E10)  are  the  maximum 
rms  differences  (percent)  allowed  between  successive  iteration 
of  the  retrieved  gas  temperature  and  concentration,  respective- 
ly, for  convergence  to  be  deemed  complete.  Both  criteria  must 
be  satisfied. 

4.  Plume  Data  Card.  The  card  name  is  PLMDAl'A.  LI  (format  E10)  is 
the  distance  (cm)  from  the  nozzle  exit  plane  to  the  observation 
scanning  plane.  L2  (format  E1U)  is  the  distance  (cm)  from  the 
scanning  plane  to  the  end  of  the  plume.  TN  (format  E10)  is  the 
temperature  (K)  of  the  nozzle  exit  plane  disc,  and  EN  (format 
E10)  is  its  emi8sivity.  These  data  (and  thus  the  card)  are 
required  only  if  MODE^BOTH  on  the  calculation  data  card  (card 
type  2). 

5.  Miscellaneous  Data  Card.  The  card  name  is  NDATA.  NSCAT  (format 
110)  is  the  number  of  equal-length  segments  that  a scattering 
line  of  sight  is  divided  into  for  numerical  integration.  Its 
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maximum  allowed  value  is  100.  NAZ  (format  110)  Is  the  number  of 
Intervals  over  which  the  360°  azimuthal  angle  Integration  Is 
performed  in  integrations  over  solid  angle.  There  is  no  dimen- 


sional limit  to  NAZ.  These  data  (and  thus  the  card)  are 
required  only  in  M0DE=*B0TH  on  the  calculation  data  card  (card 
type  2). 

Iteration  Listing  Card.  The  card  name  is  LIST.  Normally,  the 
intermediate  results  for  current  gas  temperature  and  concentra- 
tion, and  transverse  profiles  generated  at  each  iteration  are 
not  listed.  The  LIST  control  card  can  be  used  to  obtain  them. 
Listing  Suppression  Card.  The  card  name  is  NOLIST.  Listing  of 
intermediate  iteration  results  can  be  suppressed  with  the  NOLIST 
card  if  listing  was  enabled  in  previous  runs  with  a LIST  card 
(card  type  6). 

Gas  Band  Model  Parameters  Card.  The  card  name  is  GPARAM.  This 
card  calls  for  the  read  in  of  band  model  parameters  for  the  gas 
species  of  interest.  The  required  deck  structure  is  shown  in 
Fig.  14. 

Transverse  Data  Card.  The  card  name  is  ZDATA.  This  card  calls 
for  the  read  in  of  transverse  data  from  the  transverse  data 
file,  which  must  be  attached  to  the  GASIC  job  run  as  TAPE3.  The 
total  radiance  N(z),  total  transmittance  t(z),  particle-only 
radiance  Np(z),  and  particle-only  transmittance  t^(z)  functions 
are  r«'ad  In,  but  the  latter  two  functions  are  used  only  if 
MODE*FOOB  on  the  calculation  data  card  (card  type  2).  Also, 
although  it  is  not  a transverse  profile,  but  a radial  one,  the 


total  gas  pressure  profile  p (r)  is  read  in  at  this  point  from 

© 

the  transverse  data  file. 

10.  Radial  Data  Card.  The  card  name  is  RDATA.  This  card  calls  for 
the  read  in  of  radial  data  required  in  the  M0DE*BQTH  mode  of 
inversion.  If  MODE* BOTH,  the  radial  data  (and,  hence,  the  KDATA 
card)  are  not  required.  The  radial  data  required  by  the 
M0DE»B0TH  inversion  are  the  volume  absorption  and  scattering 
cross  sections  a(r)  and  8(r),  respectively,  the  scattering  phase 
function  p(r,0),  and  the  particle  temperature  Tp(r).  If  AFILii 
(format  110)  has  the  value  2,  the  a(r)  data  is  obtained  from 
the  a_  priori  radial  data  file  TAPE2.  If  AFILE  has  the  value 
4,  a(r)  is  obtained  from  the  retrieved  radial  data  file  TAPE4. 
BFILE,  PFILE,  and  TPFILE  control,  respectively,  the  input  of 

3(r),  p(r,8),  and  Tp(r)  in  the  same  way.  The  appropriate  files 
must  be  attached  to  the  job. 

11.  Execution  Card.  The  card  name  is  RUN.  When  this  card  is  en- 
countered, computations  are  begun  using  the  data  entered  up  to 
that  point,  an  output  listing  of  results  is  made,  and  the  re- 
trieved functions  Tg(r)  and  Cg(r)  are  written  to  the  retrieved 
radial  data  file  TAPE4. 
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4.  EXAMPLE  RUNS 


The  example  runs  illustrated  here  are  for  the  minimum  smoke  propellant 
(MSP)  plume  model  analysis  presented  in  Refs.  1 and  5.  Example  runs  are  given 
for  all  six  of  the  basic  codes.  Synthetic  data  were  produced  with  the  DPREP2 
and  EAPROF  codes,  and  immediate  retrieval  made  with  the  PARIC1,  PARIC2, 
PAKIC3,  AND  CASIO  codes.  Upon  completion  of  the  UPREP2  run,  the  a priori 
radial  data  file  TAPE2  was  made  permanent,  and  upon  completion  of  the  EAPROF 
run,  the  transverse  data  file  TAPE3  was  made  permanent.  These  files  were  then 
attached  to  the  inversion  runs  as  necessary.  All  of  the  inversions  used 
radial  data  from  the  a priori  radial  data  file  TAPE2,  and  thus  the  retrieved 
radial  data  file  was  never  made  permanent. 

The  codes  were  run  on  a CDC176  computer  with  SC0PE2.1  operating  system, 
and  compilation  was  made  with  the  FTN  compiler  operating  in  the  OPT-2  mode 
(slowest  compilation,  but  fastest  safe  execution).  For  the  examples  given 
here,  the  code  running  times  (compilation  plus  execution)  are  shown  in  Table 
6.  The  time  entered  in  the  table  for  GASIC  is  for  three  runs,  with  MODE  * 
CAS,  FOOB,  and  BOTH.  The  full  inversion  mode  (BOTH)  accounted  for  nearly  all 
of  the  time.  The  running  times  for  the  first  two  modes  were  only  ~ 3 sec 
each. 

4.1  Profile  Prediction  Runs 

The  MSP  plume  model  contains  A^Og  and  ^0,  respectively,  as  the  parti- 
cle and  gas  phases.  The  radial  gas  and  particle  temperature  profile  over  the 
IU  cm  radius  of  the  plume  Is  shown  in  Fig.  21.  The  gas  pressure,  gas  concen- 
tration, and  particle  loading  profiles  are  constant  in  radius  with  values  p * 
1 atm,  cg  - 0.15,  and  np  ■ 10^/cm^.  The  ^0  band  model  parameters  used  (Kef. 
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Table  6.  Code  Running  Times 


Code 


Time  (sec) 


DPREP2 

<1 

EAPROF 

84 

PARIC1 

<1 

PAR1C2 

2 

PARIC3 

19 

GASIC 

590 

* 


Total  for  three  Inversion  runs. 
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1)  are  shown  in  Fig.  22 


The  nonresonant,  self- broadening  parameter 


is  yq  * 0.07394  cm”Vatm.  The  efficiency  of  resonant  self-broadening  is  6.53, 
and  the  efficiency  for  foreign  gas  broadening  is  1.00.  Particle  scattering 
cross  sections  (Ref.  1)  were  computed  with  Mie  theory,  the  particle  size 
distribution  shown  in  Fig.  23,  and  index  of  refraction  m = 2.0  - 0.011.  The 
result  for  the  absorption  and  total  scattering  cross  sections  are  o^  ■ 3.20  x 

10"9  cm^  and  og  * 5.58  x 10“®  cm^.  The  differential  scattering  cross  section 
is  shown  in  fig.  24  (the  result  for  k * 0.01  is  used  here).  These  scattering 
parameters  were  assumed  to  be  constants  in  radius.  These  data  were  prepared 
on  the  ll-point  scattering  grid  8 » 0,  5,  15,  25,  35,  45,  6U,  90,  120,  150, 
and  180°.  The  manner  in  which  this  grid  covers  the  scattering  integral  inte- 
grand is  shown  in  Fig.  25.  The  laser  scattering  efficiency  function  was 

calculated  on  the  8-point  grid  8 - 0,  10,  2U,  3U,  45,  90,  135,  and  180°. 

L 

Radial  and  transverse  calculations  were  performed  by  dividing  the  plume  into 
10  equal-size  zones. 

The  input  data  for  DPk£P2  that  reflects  all  of  these  conditions  is 
listed  in  Fig.  26.  The  output  of  DPREP2  is  listed  in  Fig.  27.  These  results 
(as  saved  on  TAPE2)  are  the  principle  input  for  EAPR0F,  Listings  of  the 
actual  input  date  to  DPREP2  have  been  suppressed  by  deleting  the  PRINT  vari- 
able on  the  control  cards. 

Additional  input  data  for  EAPR0F  are  listed  in  Fig.  28,  Calculations 
ware  performed  for  gas-only,  particle-only,  and  coupled  gas-plus-particle 
conditions.  The  Lorentz  lineshape  and  Curtis-Godson  nonuniformity  approxima- 
tion were  used  in  the  gas  band  model.  Scattering  lines  of  sight  were  divided 
into  10  segments  for  numerical  integration.  Azimuthal  angle  integration  was 
performed  with  a 16-segment  grid.  The  distance  from  the  exit  plane  to  the 
observation  plane  was  taken  as  3 cm,  and  the  distance  from  the  observation 
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Fig.  26.  DPREP2  Input  Data  Listing  (Sheet  2 of  2) 
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Fig.  27.  DPREP2  Standard  Output  Listing  (Sheet  11  of  11) 


plane  to  the  end  of  the  plume  was  fixed  at  57  cm.  The  exit  plane  was  modeled 
as  a flat  disc  with  uniform  temperature  T * 800  K and  emlssivity  e » 0.75. 

The  output  of  EAPROF  is  listed  in  Fig.  29.  These  results  (as  saved  on 
TAPE3)  are  the  synthesized  input  transverse  data  for  the  inversion  codes. 

4.2  Retrieval  Runs 

The  input  data  for  the  PARIC1,  PARIC2,  PARIC3,  and  GASIC  runs  are 
listed  in  Figs.  30,  32,  34,  and  36,  respectively.  For  all  runs,  calculation 
condition  (e.g.,  number  of  zones,  plume  size,  lineshape)  required  a repeat  of 
the  conditions  used  in  the  prediction  run  using  DPREP2  and  EAPROF.  Any  radial 
data  required  was  always  obtained  from  the  ji  priori  radial  data  file  TAPE2 
rather  than  from  the  retrieved  radial  data  file  TAPE4.  Listing  of  input 
radial  and  transverse  data  was  suppressed.  For  PARIC2,  PARIC3,  and  GASIC,  a 
convergence  criterion  of  0.1%  within  30  iterations  was  set.  For  GASIC,  this 
criterion  applied  to  both  temperature  and  concentration.  Also  for  GASIC,  a 
multiple  run  was  made  for  each  of  the  three  retrieval  modes.  The  outputs  of 
the  four  runs  are  listed  in  Figs.  31,  33,  35,  and  37. 
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Fig.  29.  EAPROF  Standard  Output  Listing  (Sheet  2 of  13) 


xxxxxxxxxxxxxxxxxxxxLASER  SCATTERING  FUNCTIONSxxxxxxx 


CM  CM  CM  CM  CM  M KIWI 

OOOODDOO 

I I I I I I I I 
ill  m n m n 1 1 1 1 1 1 1 j 1 1 f 
U.  (Mrs.i'si'sKi—unr-i 

o rH  fo  eof^M  so  eo 
«0f*1CMi-i.-i<j-roin 


u 

Ul  oooooooo 

< ii  O oooooooo 

K w OOOOOOOO 

K X UJ  

c uu  _j  ooootnotno 

k Q O CM  rn  cn  m eo 

< Z Z Hrl 

( l-H  <t 

UJ 

01  X 

o£  ui 

uj  Q <-■  CM  ro  in  se  i>>  eo 


Fig.  29.  EAPROF  Standard  Output  Listing  (Sheet  3 of  13) 
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Fig.  29.  EAPROF  Standard  Output  Listing  (Sheet  6 of  13) 
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101 


3 unnA K Y INrUT  LISTING  kxxxxxkxxxkxx 


X 

X 

X 

X 

X 

X 

* 

* 

* 

* 

X 

X 

X 

* 

* 

* 

* 

* 

* 

X 

X 

X 

X 

* 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

if) 


rj 

in 

uu 

0£ 


\n 

o' 

LU 

> 


UJ 

—I 

u 


CxZ 

< 

ft. 

X 

X 

X 

X 

X 

X 

X 

X 

X fOrorOrOfOrOrOK>K)rOK) 

X A ooooooooooo 

X rH  I I I I I I I I I I I 

X I UJ  LU  UJ  LU  LU  LU  LU  UJ  UJ  LU  LU 

X SZ  ooooooooooo 

X O ooooooooooo 

X w (Tv  os  (7s  ffs  CT\  ^ C7>  O'  O'  ^ 

X o 

x ininifiimmAimmnimn 

x 
X 
X 

X OOOOOOOOO*— f 

X OOOOOOOOOO 

X + + + + *f  + + + *f  + 

X TZ  UJ  UJ  UJ  UJ  UJ  UJ  UJ  UJ  UJ  UJ 

X O OOOOOOOOOO 

X w OOOOOOOOOO 

X OC  OOOOOOOOOO 

X 

x o«-iCMr*>  tf*irw>r^«0<*r-4 


X 

UJ 

D rH  cj  ro  xfiri 'oj  eo  cn  o fH 
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APPENDIX 


REVISED  ABEL  INVERSION  PROCEDURE 


A transverse  function  f(z)  defined  in  terms  of  a radial  function  g(r) 
by 


<> 


f(z) 


R 


rdf 


2 / g(r)  - 

z [r  - z J 


2 _2~\  V2 


can  be  inverted  to  obtain  g(r)  with  the  Abel  inversion 


g(r) 


1 dF(r) 
2itr  dr 


(Al) 


(A2) 


where 


£ 

FCr)  - 2 / Hz)  I“fal . „ . (A3) 

(*2  - r2) 

In  the  previously  developed  gas-only  inversion  code  EMABIC  (Ref.  3),  evalua- 
tion of  the  integral  form 

• 2 J «*)  -«■— ATT  (A4) 

y (*  - y i 1 

used  in  both  Eqs.  (Al)  and  ( ii)  was  performed  with  i quadrature  formula  of  the 
form 


F - A KKSt*  E c Ck.n)  C + o <k)  G j (A5) 

n«k+l 

I 
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(with  the  summation  deleted  for  k * N and  = 0 for  k * N + 1).  In  Eq.  (A5), 
the  a's  are  quadrature  coefficients  and  and  Gn  are  function  values  at  the 

quadrature  grid  points.  R is  the  cylindrical  source  radius,  N is  the  number 
of  equal-sized  slabs  assumed  in  the  quadrature  approximation,  and  A - K/N. 
The  a coefficients  in  Eq.  A5  were  derived  with  the  assumption  that  G(x)  in  Eq. 
(A4)  varied  as 


G(x)  A + Bx^  (Ao) 

between  grid  points. 

The  use  of  Eq.  (AS)  to  evaluate  the  integrals  in  both  Eqs.  (Al)  and 
(Aj)  is  somewhat  inconsistent.  For  example,  if  the  quadratic  form  of  Eq.  (Ab) 
is  used  to  approximate  g(r)  in  the  last  radial  zone  (i.e.,  between  points  N 
and  N + 1),  then  Eq.  (Al)  predicts  that  f(z)  should  have  the  form 

f(z)  = 2(A  + Bz  ) (R  - Z J + — j (R  - z ) (A/) 

in  the  last  transverse  zone.  In  using  Eq.  (A3)  for  inversion,  however,  t(z; 
is  implicitly  assumed  to  vary  in  the  last  zone  as 

f(z)  = a + bz^,  (AB) 

[i.e.,  since  Eq.  (A5)  is  usedj.  Although  a and  b can  be  selected  so  that  Eq. 
(AB)  fits  Eq.  (A7)  reasonably  well,  it  is  not  good  enough  to  give  an  accurate 
recovery  for  g(r).  This  lack  of  ability  to  make  a good  fit  occurs  mainly  in 
the  outer  zone  and  is  caused  by  the  fact  that  the  correct  form,  Eq.  (A7),  for 
f(z)  is  very  sharply  bent  in  the  region  around  z = R.  The  derivative  t'(z) 
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given  by  Eq.  (A7)  goes  to  -*»  at  z * R while  tne  implicitly  used  form  Eq.  (AB) 
could  never  achieve  this  slope.  For  inner  zones,  the  form  Eq.  (AB)  can  be 
very  much  more  accurately  chosen  to  fit  the  true  form  of  f(z). 

The  accuracy  of  recovery  in  the  last  zone  is  further  reduced  because  a 
backward  difference  approximation  using  only  the  points  for  F(r)  at  N and 
N + 1 must  be  used  to  evaluate  g(r)  from  Eq.  (A2).  For  inner  zones,  a more 
accurate  three-point,  central-difference  formula  can  be  used* 

All  told,  the  formulation  as  presented  so  far  is  adequate  tor  profile 
generation  and  for  inversion  except  in  the  outermost  (and  sometimes  the 
second-to-last)  zone. 

This  problem  did  not  show  up  in  the  gas-only  code  EMAB1C  because  a 
straightforward  Abel  inversion  of  transverse  profiles  was  never  made.  Kather, 
Abel  inversion  was  made  on  the  difference  of  two  transverse  profiles,  which, 
in  the  limit  of  convergence,  approached  zero.  In  the  particle-only  inversion 
code  PARIC1,  however,  the  principal  output  y(r)  is  obtained  as  a direct  Abel 
Inversion  of  -in  and  serious  retrieval  errors  were  obtained. 

The  first  approach  to  solving  this  problem  was  to  retain  the  quaratic 
form  of  Eq.  (Ab)  for  g(r)  in  evaluating  f(z)  with  Eq.  (Al)  [i.e.,  in  evalu- 
ating with  Eq.  ( A5)J  but  to  derive  new  a-coef f icients  based  on  Lhe  varia- 
tion of  Eq.  ( A7 ) for  evaluating  F(r)  in  Eq.  (A3).  This  approach  was  soon 
abandoned  because  the  resulting  integrals  were  too  complicated. 

The  final  approach  was  a straightforward  back-solution  of  Eq.  (Ab). 
Assume  that  Gn  (n  * k + 1,...,  N + 1)  have  been  found.  Then,  from  Kq.  (ab), 
(k  = i,...,N)  can  be  obtained  by 


G 


k 


1 

o^k) 


I a (k,n)  G - <1  (k)  G ] 
n-k+l 


(A9) 
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(Note:  the  suaaation  term  is  deleted  for  k - N).  The  solutions  for  at  k * 
N + 1 is  determined  by  fitting  the  known  true  fora  Eq.  (A7)  to  the  points 
and  Fjj  by 


r«  - 1 - + »*%  - j)  (r2-  4 - oin  + “ <»2  - 4 - 1>3/2 


F„  - 2(A  + B.2)  (R2  - Br2)1/2  + (R2  - b2)3/2. 


(A10) 


With  Zjj  _ j * R - 2A  and  Zj|  ■ R - A,  these  two  equations  can  be  solved  for  A 
and  B.  When  determined,  they  are  substituted  into  Eq.  (A6)  and  Gjj  + j deter- 
mined by  evaluation  at  x * R.  The  result  is 


'N  + 1 


8(N  - 1)3/2  Fn  - (2N  - 1>3/2  Fn  _ j 
4A(2N  - 1)  / (2N  - 1)  (N  - 1) 


Equations  (A9)  and  (A10)  then  defined  an  inversion  procedure  for  obtaining  the 
radial  function  G(r)  from  the  transverse  function  F(z). 


